\documentclass[11pt]{article}

\usepackage{fullpage}
\usepackage{rotating}   
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{fancyhdr}
\usepackage{algorithm}
\usepackage{algorithmic}
\usepackage{bm}
\usepackage{listings}
\usepackage{graphicx}
\usepackage{caption2}
\usepackage{subfigure}
\usepackage{float}
\usepackage{extpfeil}
\usepackage{color}
\usepackage[usenames,dvipsnames]{xcolor}


\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{conjecture}[theorem]{Conjecture}
\newtheorem{remark}[subsection]{Remark}

%%
\newcommand\numberthis{\addtocounter{equation}{1}\tag{\theequation}}

%% define new symbols
\def\bx{\bm{x}}
\def\bb{\bm{b}}
\def\ba{\bm{a}}
\def\bc{\bm{c}}
\def\bf{\bm{f}}
\def\by{\bm{y}}
\def\bu{\bm{u}}
\def\bv{\bm{v}}
\def\BW{\bm{W}}
\def\BA{\bm{A}}
\def\bz{\bm{z}}
\def\BZ{\bm{Z}}
\def\BH{\bm{H}}
\def\BL{\bm{L}}
\def\BU{\bm{U}}
\def\BV{\bm{V}}
\def\BB{\bm{B}}
\def\BC{\bm{C}}
\def\BD{\bm{D}}
\def\BE{\bm{E}}
\def\BW{\bm{W}}
\def\BQ{\bm{Q}}
\def\BG{\bm{G}}
\def\BA{\bm{A}}
\def\BX{\bm{X}}
\def\BY{\bm{Y}}
\def\BQ{\bm{Q}}
\def\BI{\bm{I}}
\def\BR{\bm{R}}

%% define new brackets
\def\la{\left\langle}
\def\ra{\right\rangle}
\def\ln{\left\|}
\def\rn{\right\|}
\def\lb{\left(}
\def\rb{\right)}
\def\lsb{\left[}
\def\rsb{\right]}
\def\lcb{\left\{}
\def\rcb{\right\}}

%%
\DeclareMathOperator*{\argmin}{arg\,min}
\DeclareMathOperator*{\argmax}{arg\,max}

%%
\title{Homework XII}
\author{Name: Shao Yanjun, Number: 19307110036}


\begin{document}
\maketitle

%------------------------------------
\begin{abstract}
This is Daniel's homework of  "Numerical Algorithms with Case Studies II".
\end{abstract}
%-------------------------------------
%=====================
\section{Problems}
\paragraph{Q1}
Implement the code and examine the truncation error of RK4 on loglog() plots.
\begin{figure}[H]
	\centering
	\subfigure[truncation error]{
		\includegraphics[width=0.6\linewidth]{error.jpg}
	}
\end{figure}
Most of the truncation error satisfy $O(h^4)$ level, while further accuracy cannot be reached because of rounding error.
\paragraph{Q2}
In this implicit method, we will have the following loop invariant in each iteration.
\begin{align}
	u_{k+1}&=u_k+h(\theta f(t_k+h,u_{k+1})+(1-\theta)f(t_k,u_k))\\
	&=u_k-h\theta\lambda u_{k+1}-h(1-\theta)\lambda u_k
\end{align}
So after each iteration, we will have,
\begin{align}
	u_{k+1}=\frac{1-h(1-\theta)\lambda}{1+h\theta\lambda}u_k
\end{align}
And to make $u_k\rightarrow0$ when $k\rightarrow\infty$, we have to let,
\begin{align}
	\frac{1-h(1-\theta)\lambda}{1+h\theta\lambda}<1
\end{align}
Which gives $h>0$.\\
However, this is not sufficient for convergence, because we must take a good step $h$ to make sure every implicit equation converges within each iteration. That gives the condition that $x=(1-h(1-\theta)\lambda)y-h\theta\lambda x$ must converge.
\begin{align}
	\|x_k+1-x_k\|&=h\theta\lambda\|x_k-x_{k-1}\|\\
	h&<\frac{1}{\theta\lambda}
\end{align} 
Therefore, $h\in (0,\dfrac{1}{\theta\lambda})$. However, numerical studies show that there could be an even tighter upper bound for stepsize.
\paragraph{Q3}
The stationary concentrations are approximately $X=9.7754mg\cdot L^{-1}, C=9.7749mg\cdot L^{-1}
,S=81.4495mg\cdot L^{-1}$
\begin{figure}[H]
	\centering
	\subfigure[solution]{
		\includegraphics[width=0.6\linewidth]{res.jpg}
	}
\end{figure}
\paragraph{Q4}
\subparagraph{(a)}
We can substitute $x(t)$ with a guess of $x(t)=e^{\lambda t}$ and that gives,
\begin{align}
	(m\lambda^2+k)(e^{\lambda t})=0
\end{align}
Since $e^{\lambda t}\ne 0$, we have $\lambda=i\sqrt{\frac{k}{m}}$. So the solution will be a linear combination of $e^{i\sqrt{\frac{k}{m}}t}$ and $e^{-i\sqrt{\frac{k}{m}}t}$.
\begin{align}
	x(t)=C_1cos(\sqrt{\frac{k}{m}}t)+C_2sin(\sqrt{\frac{k}{m}}t)
\end{align}
With condition that $x(t)=x_0$, $x'(t)=v_0$, we will have,
\begin{align}
	x(t)=x_0cos(\sqrt{\frac{k}{m}}t)+v_0sin(\sqrt{\frac{k}{m}}t)
\end{align}
\subparagraph{(b)}
Let $y(t)=x'(t)$, and the second-order ordinary differential equation is transformed into a two-dimensional system of first-order coupled ordinary differential equations.
\begin{equation}
	\left\{
	\begin{aligned}
		&y'(t)=-\frac{k}{m}x(t)\\
		&x'(t)=y	
	\end{aligned}
	\right.
\end{equation}
This is easy to solve with numerical methods.
\subparagraph{(c)}
Apply four methods, and we can see that Trapezoidal and Runge-Kutta 4 are two stable methods that can be used to solve this system. Forward-Euler method seems to go too big, and Backward-Euler method seems to shrink to small.
\begin{figure}[H]
	\centering
	\subfigure[Euler method]{
		\includegraphics[width=0.4\linewidth]{with euler.jpg}
	}
	\subfigure[other methods]{
		\includegraphics[width=0.4\linewidth]{without euler.jpg}
	}
\end{figure}
And also shown below is the phase diagram of the solution.
\begin{figure}[H]
	\centering
	\subfigure[phase diagram]{
		\includegraphics[width=0.6\linewidth]{phase diagram.jpg}
	}
	
\end{figure}
Since both Euler methods are $O(h)$ methods, when $h$ is not so small, the error for $x(t)$ and $x'(t)$ accumulates very fast. What's more, as $x(t)$ and $x'(t)$ are dependent on each other, the final solution slips away from accurate ones quickly.
%-------------------------------------
%=====================
\end{document}
